#!/bin/bash -e

function info() {
echo Usage: `basename $0` '[-l] bam chr:pos'
exit 1
}

while getopts  ":p:f:i:l:" opt; do
    case  $opt  in
        p) out_prefix=$OPTARG;;
        f) suffix=$OPTARG;;
        i) interval=$OPTARG;;
        l) len=$OPTARG;;
        *) info;;
    esac
done
shift $(($OPTIND - 1))


if [ $# -lt 2 ]; then info; fi

. $var

len=${len:=1}
chr_pos=`echo $2|sed 's/\:/_/'`

samtools tview -dT -p$2 $1 $ref_genome > $out_prefix.tview.$chr_pos.$len

awk -v len=$len '{s=substr($0, 1, len); print s}' $out_prefix.tview.$chr_pos.$len > $out_prefix.$chr_pos.$len

head -2 $out_prefix.$chr_pos.$len |tail -1 > $out_prefix.$chr_pos.$len.ref
sed '1,2d; /^  *$/d' $out_prefix.$chr_pos.$len|sort|uniq -c|sed 's/^  *//; s/ /\t/' |sort -k1nr > $out_prefix.$chr_pos.$len.stat

